package fem2;

import inf.jlinalg.Array1DMatrix;
import inf.jlinalg.ArrayVector;
import inf.jlinalg.BLAM;
import inf.jlinalg.IMatrix;
import inf.jlinalg.IVector;

import java.util.ArrayList;

import math2.FunctionBasisOnRnToR;
import math2.LinAlgUtilities;
import math2.PolynomialBasisOnRnToR;
import fem2.enu.MeshPartType;

/**
 * represents a patch in the mesh. It is used for stress recovery
 * 
 * @author hbui
 * 
 */
public class Patch {

	private ArrayList<MeshPart> mpList;
	private int p = 0; // order of the mesh part in the patch
	private boolean bilinear = false;

	public Patch(ArrayList<MeshPart> mpList) {
		this.mpList = mpList;

		for (int i = 0; i < mpList.size(); i++) {
			FunctionBasisOnRnToR basis = mpList.get(i).getBasis();
			int p;
			if (basis instanceof PolynomialBasisOnRnToR) {
				p = ((PolynomialBasisOnRnToR) basis).getP();
			} else {
				/*
				 * TODO: check this
				 */
				p = 5; // set a large p to account for accuracy
			}
			if (this.p == 0) {
				this.p = p;
				if (mpList.get(i).getType() == MeshPartType.QUAD4
						|| mpList.get(i).getType() == MeshPartType.HEX8) {
					bilinear = true;
				} else {
					if (bilinear) {
						throw new Error("mesh part was not the same bilinear type. It's invalid");
					}
				}
			} else if (p != this.p) {
				throw new Error("mesh part in the same patch must be of the same order");
			}
		}
	}

	public int countMeshPart() {
		return mpList.size();
	}

	public MeshPart getMeshPart(int i) {
		return mpList.get(i);
	}

	/**
	 * get sampling depth for left hand side matrix to ensure the number of
	 * sampling points is more than number of terms in interpolation matrix
	 */
	private int getSamplingDepthForA() {
		int dim = mpList.get(0).getDimension();
		if (this.p == 1) {
			if (mpList.size() < dim + 1) {
				return this.p + 1;
			} else {
				return this.p;
			}
		}
		return this.p;
	}

	/**
	 * compute interpolation vector P
	 * 
	 * @param x
	 * @return
	 */
	public IVector computeInterpolationVector(double... x) {
		return computeInterpolationVector(p, bilinear, x);
	}

	/**
	 * compute interpolation vector P
	 * 
	 * @param p
	 * @param bilinear
	 * @param x
	 * @return
	 */
	public static IVector computeInterpolationVector(int p, boolean bilinear, double... x) {
		int dim = x.length;
		if (p == 1) {
			if (dim == 2) {
				return new ArrayVector(1.0, x[0], x[1]);
			} else if (dim == 3) {
				return new ArrayVector(1.0, x[0], x[1], x[2]);
			} else {
				throw new Error("invalid dimension");
			}
		} else if (p == 2) {
			if (bilinear) {
				if (dim == 2) {
					return new ArrayVector(1.0, x[0], x[1], x[0] * x[1]);
				} else if (dim == 3) {
					return new ArrayVector(1.0, x[0], x[1], x[2], x[0] * x[1], x[1] * x[2], x[0]
							* x[2]);
				} else {
					throw new Error("invalid dimension");
				}
			} else {
				if (dim == 2) {
					return new ArrayVector(1.0, x[0], x[1], x[0] * x[0], x[0] * x[1], x[1] * x[1]);
				} else if (dim == 3) {
					return new ArrayVector(1.0, x[0], x[1], x[2], x[0] * x[0], x[1] * x[1], x[2]
							* x[2], x[0] * x[1], x[1] * x[2], x[0] * x[2]);
				} else {
					throw new Error("invalid dimension");
				}
			}
		} else {
			throw new Error("interpolation function with higher order than 2 is not supported yet");
		}
	}

	/**
	 * compute left hand side matrix A
	 * 
	 * @param A
	 * @return
	 */
	private IMatrix computeA(IMatrix A) {
		IntegrationManager integrator = IntegrationManager.getInstance();
		int d = getSamplingDepthForA();
		for (int i = 0; i < mpList.size(); i++) {
			MeshPart mp = mpList.get(i);
			MeshPartType mpType = mp.getType();
			/*
			 * integrate over an element in a patch
			 */
			int ng = integrator.getNumIntegrationPoints(mpType, d);
			for (int j = 0; j < ng; j++) {
				double[] xi = integrator.getIntegrationPoint(mpType, d, j);
				double w = integrator.getIntegrationWeight(mpType, d, j);
				double[] X = mp.XAt(xi);
				IVector P = computeInterpolationVector(X);
				// Debugger.watch("P = ", P);
				BLAM.multiply(w, P, P, A);
			}
		}
		return A;
	}

	/**
	 * solve for coefficients
	 * 
	 * @param values
	 *            the values at respective Gauss point (must be in the
	 *            appropriate order)
	 * @return
	 */
	public IVector solveForCoefficients(double... values) {
		int n = 0;
		IVector b = null;
		/*
		 * compute right hand side
		 */
		int cnt = 0;
		for (int i = 0; i < mpList.size(); i++) {
			MeshPart mp = mpList.get(i);
			/*
			 * integrate over an element in a patch. To integrate b, the number
			 * of sampling points is get equal to number of integration points
			 * on the element
			 */
			int ng = mp.getElement().getNumIntegrationPoints();
			for (int j = 0; j < ng; j++) {
				double[] xi = mp.getElement().getIntegrationPoint(j);
				double w = mp.getElement().getIntegrationWeight(j);
				double[] X = mp.XAt(xi);
				if (b == null) {
					b = computeInterpolationVector(X);
					n = b.getSize();
					BLAM.scale(w * values[cnt++], b);
				} else {
					IVector P = computeInterpolationVector(X);
					BLAM.add(w * values[cnt++], P, b);
				}
			}
		}
		/*
		 * compute left hand side
		 */
		IMatrix A = Array1DMatrix.createColumnMajor(n, n);
		A = computeA(A);
		// Debugger.watch("A = ", A);
		/*
		 * solve for coefficients
		 */
		LinAlgUtilities.backSolve((Array1DMatrix) A, ((ArrayVector) b).getValues());
		return b;
	}

	@Override
	public boolean equals(Object o) {
		if (this == o) {
			return true;
		}
		if (!(o instanceof Patch)) {
			return false;
		}
		Patch p = (Patch) o;
		if (p.countMeshPart() != mpList.size()) {
			return false;
		} else {
			for (int i = 0; i < p.countMeshPart(); i++) {
				if (!mpList.contains(p.getMeshPart(i))) {
					return false;
				}
			}
			return true;
		}
	}
}
